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Abstract  -  We  describe  a  primal-dual  interior  point  algorithm  for  convex  qua¬ 
dratic  programming  problems  which  requires  a  total  of  0(a^L)  arithmetic  opera¬ 
tions.  Each  iteration  updates  a  penalty  parameter  and  finds  an  approximate 
Newton’s  direction  associated  with  the  Kuhn-Tucker  system  of  equations  which 
characterizes  a  solution  of  the  logarithm  barrier  function  problem.  This  direction 
is  then  used  to  find  the  next  iterate.  The  algorithm  is  based  on  the  path  follow¬ 
ing  idea.  The  total  number  of  iterations  is  shown  to  be  of  the  order  of  O(^fnL). 

Key  Words .  -  Interior-point  methods;  Convex  Quadratic  Programming; 
Karmarkar’s  algorithm.  Polynomial-time  algorithms,  Barrier  function.  Path  fol¬ 
lowing. 


1.  Introduction 

In  Monteiro  and  Adler  [12],  an  algorithm  to  solve  Linear  Programming  problems  has 
been  presented  which  converges  in  0(4nL)  iterations  with  an  average  number  of  0(n 2S) 
arithmetic  operations  per  iteration.  In  the  last  section  of  that  paper,  the  authors  observed 
that  the  same  techniques  could  be  extended  to  solve  convex  Quadratic  Programming  prob¬ 
lems  in  at  most  0(\^nL)  iterations.  The  objective  of  the  present  paper  is  to  present  the 
details  of  the  algorithm  mentioned  in  [12]  as  applied  to  convex  Quadratic  Programming 
problems. 

Quadratic  Programming  (  QP  )  problems  share  many  of  the  combinatorial  properties  of 
Linear  Programming  (  LP  )  problems.  Based  on  these  properties,  algorithms  extending  the 
simplex  method  have  been  devised  to  solve  QP  problems.  However,  in  the  worst  case,  these 
algorithms  may  converge  in  an  exponential  number  of  steps. 


Polynomial-time  algorithms  for  convex  Quadratic  Programming  problems  based  on  the 
ellipsoid  method  were  presented  in  [1]  and  f  10). 

Recently,  with  the  advent  of  the  new  interior  point  algorithm  by  Karmarkar  (7]  for 
solving  LP  problems,  some  attention  has  been  devoted  to  study  classes  of  problems  that  can 
be  solved  by  interior  point  algorithms  in  polynomial  time.  Ye  and  Tse  [15]  present  an  inte¬ 
rior  point  algorithm  for  solving  convex  QP  problems  based  on  Karmarkar’s  projective 
transformation.  Their  algorithm  is  shown  to  converge  in  at  most  O(nL)  iterations  with  a 
computational  effort  of  0(n3L2)  arithmetic  operations  per  iteration.  Thus  overall  their 
algorithm  involves  0(n4L3)  arithmetic  operations  in  the  worst  case. 

The  algorithm  discussed  in  this  paper  is  based  on  the  logarithm  barrier  function 
method  and  on  the  idea  of  following  the  path  of  minimizers  for  the  logarithm  barrier  family 
problems,  that  is,  the  so  called  "central  path".  This  path  has  been  extensively  studied  in 
Bayer  &  Lagarias  [2]  and  Meggido  [11].  The  logarithm  barrier  function  approach  is  usually 
attributed  to  Frisch  [4]  and  is  formally  studied  in  [3]  in  the  context  of  nonlinear  optimiza¬ 
tion.  Algorithms  for  LP  problems  based  on  following  the  central  path  have  been  presented 
in  [13],  [14],  [6]  and  [12].  The  breakthrough  in  this  line  of  research  was  obtained  by  Rene- 
gar  [13],  who  was  the  first  to  achieve  a  speed  of  convergence  of  Oi^fnL)  iterations,  where 
each  iteration  involves  0{ny)  arithmetic  operations.  His  algorithm  is  based  on  the  method 
of  centers  following  the  central  path.  Subsequently,  Vaidya  [14]  improved  Renegar’s  com¬ 
plexity  to  a  total  of  0(n3L)  arithmetic  operations  using  the  same  approach  of  the  method  of 
centers.  His  algorithm  converges  with  the  same  order  of  iterations  as  in  Renegar’s  algo¬ 
rithm,  however,  he  showed  that  the  average  number  of  arithmetic  operations  per  iteration 
can  be  bounded  by  0(n 2  *).  Independently,  an  equivalent  complexity  was  also  obtained  by 
Gonzaga  [6],  using  the  logarithm  barrier  function  approach.  Both  Vaidya’s  and  Gonzaga’s 
algorithms  arc  primal  algorithms.  Kojima  ct  al.  [8]  presented  a  primal-dual  algorithm  based 
on  the  logarithm  barrier  function  method  and  the  primal-dual  framework  described  in  [11], 
Their  algorithm  is  shown  to  converge  in  al  most  O(nL)  iterations  with  a  computational 


effort  of  0(n3)  arithmetic  operations  per  iteration,  resulting  in  a  total  of  0(n*L)  arithmetic 
operations.  Based  on  Kojima  et  al.  [8]  and  on  Gonzaga  [6],  a  primal-dual  algorithm  con¬ 
verging  in  Oi'fnL)  iterations,  in  the  worst  case,  with  an  average  computational  effort  per 
iteration  of  0(n **)  arithmetic  operations  was  presented  by  Monleiro  and  Adler  [12].  As 
mentioned  above,  the  current  paper  is  an  extension  of  these  thechniques  as  applied  to  con¬ 
vex  Quadratic  Programming  problems  and  it  achieves  a  complexity  similar  to  the  Linear 
Programming  case  presented  in  [12]. 

Our  paper  is  organized  as  follows.  In  section  2,  we  present  some  theoretical  back¬ 
ground.  In  section  3,  we  present  the  algorithm.  In  section  4,  we  prove  results  related  to  the 
convergence  properties  of  the  algorithm  and  we  also  describe  the  updating  scheme  that 
leads  to  a  reduction  in  the  average  number  of  arithmetic  operations  per  step.  In  section  5, 
we  discuss  how  to  initialize  the  algorithm.  In  section  6,  we  conclude  the  paper  with  some 
remarks. 

2.  Theoretical  Background 

In  this  section,  we  briefly  review  some  theoretical  results  pertinent  to  the  present 
work.  A  detailed  discussion  of  these  results  can  be  found  in  [11],  We  consider  the  convex 
quadratic  programming  problem  as  follows.  Let 

(Pi  min  cTx+  *  xT Q.x 

St  A\  =  /> 

x  >  0 

where  t  ,  x  are  /i-vectors,  h  is  an  m-veclor,  A  is  an  mxn  matrix  and  Q  is  a  positive  semi- 
definite  nxn  matrix.  As  for  linear  programming  problems,  the  following  fact  is  true  for 
convex  quadratic  programming  problems. 


Proposition  2.1:  If  problem  (P)  does  not  have  an  optimal  solution  then  it  must  be  either 
unbounded  or  infeasible. 

The  LagTangian  dual  problem  corresponding  to  problem  (P)  is  another  quadratic  pro¬ 
gramming  problem  given  by 

( D )  max  -  ^  vrQv  +  bTy 

s.t.  -  Qv  +  At y  +  :  =  c 

Z  >  0 

where  v  and  z  are  n-vectors  and  y  is  an  m-vector.  The  relationship  between  problems  (P) 
and  (D)  is  provided  by  the  following  result  known  as  the  duality  theorem  for  convex  qua¬ 
dratic  programming. 

Proposition  2.2:  The  follow  ing  statements  arc  true 

(a)  If  problem  (P>  is  unbounded  then  problem  ( D )  is  infeasible.  If  problem  ( D )  is 
unbounded  then  problem  (Pi  is  infeasible. 

(bi  If  problem  (Pi  has  an  optimal  solution  x“  then  there  exist  y°  and  z°  such  that  the  point 
(v.  y,  z)  =  ( x“,  y‘\  z°)  is  an  optimal  solution  of  problem  (D).  Conversely,  if  problem  (D) 
has  an  optima l  solution  then  problem  (P)  has  an  optimal  solution.  Moreover,  the  optimal 
values  of  both  problems  are  identical. 

The  complementary  slackness  version  for  convex  quadratic  programming  problems  is 
as  follows. 

Proposition  2.3:  If  x°  and  (v",  y‘\  z")  ate  optimal  solutions  for  problems  (P)  and  (D) 
respectively  then 

(x°)Tz°  =  0  (2.1) 

Conversely,  if  (v,  y,  z)  =  (x",  y°,  z°)  is  a  feasible  solution  of  (D)  such  that  x°  feasible  for 
(P)  and  such  that  relation  (2.1)  holds,  then  x"  and  (.r",  y",  r")  arc  optimal  solutions  of 


problems  (P)  and  ( D )  respectively. 


Tbe  algorithm  we  consider  in  this  paper  has  its  motivation  on  the  application  of  the 
logarithm  barrier  function  technique  to  problem  (P).  The  logarithm  barrier  function  method 
consists  of  a  consideration  of  the  family  of  problems 

1  n 

(P^)  min  cT x  +  xTQx  -  p  £  In 
—  .-1 


where  p  >  0  is  the  barrier  penalty  parameter.  This  technique  is  well-known  in  the  context 
of  general  constraint  optimization  problems.  One  solves  the  problem  penalized  by  the  loga¬ 
rithm  barrier  function  term  for  several  values  of  the  parameter  p,  with  p  decreasing  to  zero, 
and  the  result  is  a  sequence  of  feasible  points  converging  to  a  "solution"  of  the  original 
problem.  This  method  is  usually  attributed  to  Frisch  |4].  The  interested  reader  can  refer  to 
Fiacco  &  McCormick  [3]  for  a  detailed  discussion  of  this  technique  in  the  context  of  non¬ 
linear  constrained  optimization.  Recently  this  method  was  first  reconsidered  in  [5]  where  a 
similarity  with  Karmarkar's  algorithm  is  discussed.  Mcggido  [11]  gives  a  comprehensive 
analysis  of  the  logarithm  barrier  function  approach  as  applied  to  Linear  Programming  and 
Linear  Complementary  problems  with  positive  semi-definite  matrices. 

Before  we  can  apply  the  logarithm  barrier  function  method,  some  assumptions  on  the 
problems  (P)  and  (D)  are  necessary.  We  impose  the  following  assumptions: 

Assumption  2.4: 

(a)  The  set  S  s  j  x  e  R"  ;  Ax  =  h,  x  >  0  J  is  non  empty. 

(b)  The  set  T  =  j  (v,  v,  r)e  R”xRmxR''  ;  -  Q\  +  A  7  v  +  :  —  c  ,  r  >  0  j  is  non-empty . 


We  say  that  points  in  the  sets  5  and  T  are  interior  feasible  solutions  of  problems  (P)  and 
(D)  respectively.  TT  need  for  (a)  is  evident  since  the  logarithm  barrier  function  method 
always  works  in  the  interior  of  the  set  defined  by  the  inequality  constraints.  Assumptions 

(b)  and  (c)  arc  also  necessary  as  w  ill  become  clear  from  the  discussion  that  follows.  In  sec¬ 

tion  5,  we  will  show  how  one  can  transform  any  given  problem  to  one  satisfying  assump¬ 
tion  2.4. 

Throughout  this  paper,  we  use  the  following  notation.  If  x  =  xn)T  is  an  n- 

vector,  then  the  corresponding  capital  letter  X  denotes  the  diagonal  matrix  diag(x j,...,  xn). 
Observe  that  the  objective  function  of  problem  ( P is  a  strictly  convex  function.  This 

implies  that  the  problem  (P^)  has  at  most  one  global  minimum,  and  that  this  global 

minimum,  if  it  exists,  is  completely  characterized  by  the  Karush-Kuhn-Tucker  stationary 
condition: 

c  +  Q\  -  -  A7 \  =  0 

Ax  =  b  ,  .v  >  0 

where  e  denotes  the  n-vector  of  ones  and  y  is  the  Lagrangian  multiplier  associated  with  the 
equality  constraints  of  problem  (Pp).  By  introducing  the  n-vector  r,  this  system  can  be 
rewritten  in  an  equivalent  way  as 

(/ )  ZXe  -  pr  =  0  (2.2) 

(</ )  Ax  =  h  ,  x  >  0 

(iu)  -  Qi  +  A7  y  +  :  =  c 

A  necessary  and  sufficient  condition  for  the  problem  (PM)  to  have  a  solution  for  all  /j  >  0  is 
given  by  the  following  result. 


I 


-  7  - 


Proposition  2.5  :  Assume  ( a )  of  Assumption  2  4  holds  and  let  p  >  0  be  given  Then  prob¬ 
lem  ( P p)  has  an  optimal  solution  if,  and  only  if.  the  set  of  optimal  solutions  of  problem  ( P ) 
is  non-empty  and  bounded. 

From  this  result,  we  immediately  conclude  that  if  ( P  )  has  a  solution  for  some  p  >  0 
then  it  has  a  solution  for  all  p  >  0.  The  role  played  by  assumption  (b)  is  now  provided  by 
the  following  result. 

Proposition  2.6  :  Assume  that  problem  (P)  is  feasible.  Then  the  set  of  optimal  solutions  of 
problem  (P)  is  non-empty  and  bounded  if,  and  only  if,  assumption  (b)  holds,  that  is,  the  set 
of  interior  feasible  solutions  of  the  dual  problem  (D)  is  non-empty. 

As  a  consequence  of  the  previous  two  propositions,  we  have  the  the  following  corol¬ 
lary. 

Corollary  2.7  :  Under  assumptions  (at  and  ( b ),  problem  ( P (and  consequently  the  system 
(2.2)  )  has  a  unique  solution  x(p),for  all  p  >  0 

The  Kuhn-Tucker  system  (2.2)  provides  important  information  which  we  now  point 
out.  Assume  that  p  >  0  is  fixed  in  the  system  (2.2).  Since  x  >  0,  the  first  equation  in  the 
system  (2.2)  implies  that  r  >  0.  The  third  equation  m  (2.2)  then  implies  that  the  triple 
(a,  y,  z)  is  an  interior  feasible  solution  for  the  dual  problem  (D).  From  assumption  (c),  it 
follows  that  there  is  a  unique  y  satisfying  (2.2).  We  denote  the  unique  triple  that  satisfies 
(2.2)  by  w(p)  =  (a ip),  yip),  zip)).  Consider  the  set  H  defined  by 

W  =  |  (a,  y,  r)e  R^xR^xR"  ;  Av  =  b  .  a  >  0  ,  -Qx  +  ATy  +  :  =  c  ,  :  >  0 

Observe  that  W  is  the  set  consisting  of  the  interior  feasible  solutions  ix,  y,  z)  of  problem 
(D)  such  that  a  is  an  interior  feasible  solution  of  problem  (P).  Obviously  w (p )  is  in  the  set 
W.  The  duality  gap  at  point  we  W  is  by  definition  given  by 

g(w)  =  cta  +  a tQx  -  bT y 


which  is  simply  the  value  of  the  objective  function  of  the  problem  (P)  at  x  minus  the  value 
of  the  objective  function  of  the  problem  (D)  at  (x,  y.  ;).  From  the  definition  of  the  set  VV\ 
one  can  easily  verify  that  if  we  H’  then 

*(»••)  =  *T:  (2.3) 

In  particular,  using  the  first  equation  in  (2.2),  it  follows  that 
£(w(p))  =  np 

for  all  p  and  therefore  j? ( ))  converges  to  zero  as  p  approaches  zero.  This  implies  that 
the  objective  function  value  of  problem  (P)  at  xip)  and  the  objective  function  value  of 
problem  (D)  at  (x(p),  y (/j ).  z(p))  converge  to  the  common  optimal  value  of  problems  (P) 
and  (D).  In  fact,  the  following  stronger  result  holds  true  (  c.f  [11]  ). 

Proposition  2.8  :  Under  assumptions  (a),  (h)  and  (c).  as  /j  — »  0 ,  x(p)  and 

ixip).  yip),  zip))  converge  to  optima I  solutions  of  problems  (P)  and  (D)  respectively. 

The  following  notation  will  be  useful  later.  Let  w  =  (a,  y,  r)e  W.  We  denote  by 
/(w)  =  (/^m  ),  ...  ,/„(m  ))'  €  R"  the  n-vector  defined  by 

/,(*')  =  x,z,  ,t  =  l . n 

We  denote  by  T  the  set  (or  path)  of  solutions  w(q)  ,  p  >  0  for  the  system  (2.2),  i.e., 

T  =  |  vi(q)  =  (xip).  yip),  zip))  ;  p  >0 

The  algorithm  which  will  be  presented  in  the  next  section  is  based  on  the  idea  of  fol¬ 
lowing  this  path  F  closely.  The  path  T  will  serve  as  a  criterion  to  guide  the  points  generated 
by  the  algorithm. 


3.  The  Algorithm 


The  algorithm  presented  in  this  scumii  parallels  the  one  presented  in  [  1  2 J  lor  Linear 
Programming  problems  We  refer  the  reader  to  { H |  and  1 1 2 J  for  a  motivation  of  the  direc  - 
lions  generated  by  the  algorithm  that  we  now  describe.  The  directions  generated  by  the 
algorithm  are  determined  as  follows  Given  a  point  w  =  (»,  v,  r)  in  the  set  vL,  we  consider 
the  direction  Aw  =  (A».A\.A:,c  F'xK'  xR’'  determined  by  the  following  system  of 
linear  equations 

Z  At  +  A  A:  —  XZc  —  pc  (3.1  a) 

.4  A  i=0  (3.1  .b) 

—  {J  A  i  +  A  7  Ay  +  A;  =  ()  (3.1.C) 


where  .7  and  z  are  vectors  with  all  components  positive  and  /i  >  0  is  some  prespecified 
penalty  parameter.  When  x  =  7  and  :  =  r ,  the  direction  An  is  exactly  the  Newton's  direc¬ 
tion  associated  with  the  system  (2.2).  In  the  algorithm  given  below,  we  let  7  and  :  be 
approximations  of  the  points  x  and  :.  The  criterion  of  approximation  is  given  in  step  2  of 
algorithm  3.1.  After  some  algebra,  one  obtains  the  following  expressions  for  Ax  and  Av. 

r  /  \  -l  -i 


A  >  -  (Z  +  XQf'  I  -  A' A  7  i  A  (Z  +  Aer^'A7 


A  (Z  +  XQ)~] 


(XZe  -  fic) 


Ay  =  - 


u 


A  ( Z  4- 


A Q  f'\A7 


-  I 

A  (Z  + 


(AZe  -  /re) 


Before  describing  the  algorithm,  some  notations  are  necessary.  Let  s  denote  the  pair  of 
approximations  (7.  r  ).  We  denote  the  direction  Aw  =  (Av,  Ay,  A:)  determined  by  the  sys¬ 
tem  (3.1)  as 


Aw  ( w  ,  7,  /} ) 

in  order  to  indicate  its  dependence  on  the  point  w  =  (  v,  y,  ;),  on  the  approximation  s  and 


on  the  penalty  parameter  /i 


We  are  now  ready  lo  describe  the  algorithm.  At  the  beginning  of  the  algorithm,  we 


assume  that  an  initial  point  w°  =  (x°.  y°,  :“)  e  W  is  available  such  that  the  following  cri¬ 
terion  of  closeness  with  respect  to  the  path  T  is  satislied: 


/(h°)  -  n°e  ||  <  6fi° 


where  ||  .  ||  denotes  the  Euclidean  norm,  p“  is  a  positive  constant  and  6  =  0.1. 
We  now  state  the  algorithm. 


Algorithm  3.1  : 


Step  0)  Let  e.  W  and  p°  >  0  satisfy  (3.2).  Let  e  be  a  given  tolerance  for  the  duality 
gap.  Let 


S  :=  0.1 


y  :=  0.1 


Set  k:=0. 


Step  1)  If  g(wk)  =  xkT  :k  <  e,  stop 


Step  2)  Choose  s  =  (x ,:)  in  R"  x  R".  satisfying: 


K  -  .  . 

- <  y  ,  i  =  1 . n 

I  .v,  I 


y  ,  i=l . n 


Step  3)  Set  pt+]  :=  p*(l  -  8  l\n). 


Calculate  Avr  =  Aw(w*,  s,p  ). 
Step  4)  Set  wk*]  :=  w*  -  Aw*  . 


'Wl 
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Set  k  .  =  k  +  1  and  go  to  step  1 . 

In  the  following  sections,  we  prove  that  all  points  generated  by  algorithm  3.1  are  in 
the  set  H’  and  that  they  remain  close  to  the  path  T  in  a  sense  to  be  described  latter.  We 
also  show  that  it  terminates  in  at  most  0(4n  max(  logf'1,  logn,  logp°)  )  iterations. 
Finally,  we  present  a  suitable  choice  for  the  approximation  point  s  =  (x,  T)  (see  step  2  of 
the  algorithm  3.1)  that  will  enable  us  to  show  that  algorithm  3.1  performs  no  more  than 
0(n*  max(  logf'1,  logn.  logp")  )  arithmetic  operations  until  its  termination 


4.  Convergence  Results 

In  this  section,  we  present  convergence  results  for  the  algorithm  described  in  section 
3.  Similar  convergence  results  for  Linear  Programming  problems  are  presented  in  [12] 
Since  the  proofs  of  some  results  in  this  section  are  exactly  the  same  as  for  the  Linear  Pro¬ 
gramming  case,  the  interested  reader  is  referred  to  |J2]  for  a  detailed  discussion.  We  have 
omitted  the  proofs  of  those  results  which  follow  without  modification  from  [12). 

Let  h  =  (a,  y,  r)e  VL,  s  =  (a,  r  )e  R"x  R"  and  p  >  0.  Let  Ah  =  (Ax ,  Ay.  A:)  be 
the  direction  Aw(w,  ?,  / 5).  Consider  the  point  defined  by  w  =  w  -  Aw  .  The  next  result 
provides  expressions  for  the  product  of  complementary’  variables  /,(w  )  ,  i  =  1 ,  ...  , n. 


I 


1 »V 

m 

Pr 

r.m 

r , 

'  ^ 

!  A 

i 


'i 

s 


Proposition  4.1:  Let  w\  s  and  w  be  as  above  Then  the  following  expressions  hold 

/,(**)  =  P  +  A v, Ar,  +  (a,  -  v,)A:,  +  (:,  -  r,)Av,  (4.1 ) 


(AA)r(A;)  >  0 


(4.2) 


Proof:  Expression  (4.1)  can  be  easily  proved  using  the  definition  of  /,( w  )  and  expression 
(3. l.a).  Multiplying  expression  (3. 1  ,c)  on  the  left  by  (Ax)r,  we  obtain 


(A Ax)7" Ay  +  (Ar)TAr  -  (A x)TQA\  =  0 


(4.3) 
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Relations  (3.1.b),  (4.3)  and  the  fact  that  the  matrix  Q  is  positive  semidefinite  immediately 


imply  (4.2).  This  completes  the  proof  of  the  proposition.  □ 


We  now  state  and  prove  a  result  that  provides  bounds  necessary  to  show  that  the 


points  generated  by  algorithm  3.1  are  feasible  and  remain  close  to  the  path  T. 


Let  *  =  <  i.  >,  :)  e  W  s  -  (x.  : )  e  F^x  Fn  and  (i  >  0.  Let  Aw  =  (Ax,  Ay,  Ar) 


be  the  direction  Am*.  a  .  jj  )  We  denote  by  A/  s  (A/ , . A fn)T  the  n-vector  defined  as 


A/  =  t  Ai  ,A: ( .  ,.\\n\:n)T 


where  A  i  and  A:,  denotes  the  i"‘  coordinate  of  the  vectors  Ax  and  A z  respectively.  The 


r. e x t  result  provides  an  upper  bound  on  the  Euclidean  norm  of  the  vector  A/. 


l.emma  4.2:  Let  Ar  he  defined  c is  in  <4  4/  Then.  we  have 


.\i  n  <  ii 


where 


f  mm  s  min  ^  I  7,7,  I  ;  i  =  1,  ...  .« 


Furthermore .  have 


ii5i.-ir  <  hi 

f  min 


f  nun 


where  D  is  the  diagonal  matrix  defined  by 


D  =  (Z-'X  )Vl 


Proof:  By  equation  (3.1. a),  we  have 


a)  v:  ■ , 


m 


"IdfjfCv'j- v  j-  2* jhv’  V  "‘.'li  V'.’v'-*ks 


D''A.x  +  DA:  =  (A'Z  )~\XZ  -  fic) 


(4.10) 


From  (4.2)  it  follows  that 


(D']Ax)T(DA:)  >  0 


Using  relations  (4.10),  (4.11)  and  the  definition  of  the  Euclidean  norm,  we  obtain 
||  D-'Ax  |f  +  II  DA:  |f  <  ||  D~'Ax  |f  +  2(D^Ax)r(DA:)  +  ||  DA:  |f 


=  \\D~'Ax  +  D  Ar|| 


=  ||  (XZyv\XZe  -  fie) 


V  ~  A)~ 


<  ii/(wq  -  fie  n~ 

f  min 


Inequalities  (4.7)  and  (4.8)  follow  immediately  from  (4.12).  Also  (4.12)  implies  that 
||  D~'Ax  ||  ||  DA:  ||  <  ^  II* 

mm 

On  the  other  hand,  using  the  Cauchy-Schwar/.  inequality,  we  obtain 


Af  ||  <  £  I  Ax, A:, 

i  =  i 


=  £  I  D,; 1  A.v,  I  I  D„A:,  I 

i  =  i 


<  ||  D-]Ax  ||  ||  DA: 


(4.11) 


(4.12) 


(4.13) 


(4.14) 


Relations  (4.13)  and  (4.14)  imply  inequality  (4.5).  This  completes  the  proof  of  the  lemma. 


We  now  slate  the  key  result  to  prove  the  convergence  of  algorithm  3.1.  We  first  intro¬ 
duce  some  notation.  Given  two  vectors  u  R"  and  ie  R",  we  denote  the  Euclidean  norm 
of  the  vector  X_,(  Jf  -  x  )  by  ||  x  -  x\l  ,  i.e.. 


X  ~  *11,  =  X 

.  =  1 


(4.15) 


The  key  result  is: 


Theorem  4.3:  Let  w  =  (x,  y,  r)  e  W  and  p  >  0  satisfy 


/( w)  -  /if  ||  <  0/i 


(4.16) 


Let  f=(j?,F)e  R"  x  R"  satisfy 


\x,-x,\ 

— —  ^7.  '  =  1 . n 

I  x ,  I 


(4.17) 


I  r,  -  F,  I  ^ 

nrr~-r  • i  =  1 . " 


(4.18) 


Let  p  >  0  be  defined  as 


p  =  //(  1  —  S  l*Jn ) 


(4.19) 


Consider  the  point  w  &  (x,  y,  z)  e  R"x  Rmx  R"  defined  by 


H’ s  w  -  Aw 


(4.20) 


where  Aw  =  Aw(w,  s  ,  p).  Then  the  following  hold: 
(ai  The  point  w  is  in  the  set  W  and  satisfies 


|i  -  jt||,  <  0.28 


(4.21) 


-  rll.  <  0.28 


(4.22) 


(b)  \\f  (  h')  -  fie ||  <  6fi 


(c)  g(w)  s  xTz  <  l.li 


The  proof  of  theorem  4.3  is  exactly  the  same  as  for  the  Linear  Programming  case  and 
hence  will  not  be  given  here  (see  lemma  4.8  and  theorem  4.1  in  section  4  of  [12]  ).  As  a 
consequence  of  theorem  4.3,  we  have  the  following  corollary. 

Corollary  4.4:  All  points  wk  generated  by  algorithm  3.1  satisfy 

( a )  wk  is  in  the  set  W  .for  all  k  =  1,2....  and 

11  xt  +  '  -  xk  |L.  <  0.28 


:*  +  1  -  |L.  <  0.28 


(b)  ||  /(w1)  -  nke  ||  <  6pk  .for  all  k  =  1,2,... 

(c)  g(w*)  *  xkT:k  <  ].\npk  .  for  all  k  =  1,2,... 


where 


/  =  ai°(  1  -  8  hfn  )k  for  k  =  1.2 


Proof:  This  result  follows  trivially  by  arguing  inductively  and  using  theorem  4.3.  □ 

We  now  derive  an  upper  bound  on  the  total  number  of  iterations  performed  by  algo¬ 
rithm  3.1.  The  following  result  follows  easily  from  Corollary  4.4  and  is  proved  in  section  4 
of  [12]. 

Proposition  4.5:  The  total  number  of  iterations  performed  by  algorithm  3.1  is  no  greater 
than  k’  s  [  log(l.l/if-1  p")  Vn  /  6  ]  where  e  >  0  denotes  the  tolerance  for  the  duality 
gap  and  p°  is  the  initial  penalty  parameter. 

With  respect  to  the  data  of  problem  (P),  define 


£ 


L  =  log-)  l‘ar8cst  absolute  value  of  the  determinant) 
l  of  any  square  submatrix  of  M  j 

+  log:  (  max  I  t  j  I  )  +  log,  (  max  I  b,  I  )  +  log-,  (  m  +  n  ) 

j  i 


(4.23) 


where  M  is  the  matrix  given  by 


Q  a 

At  0 


(4.24) 


It  is  straightforward  to  verify  that  the  constant  L  is  is  less  than  two  times  the  number  of 
bus  necessary  to  represent  the  data  of  problem  (P).  The  following  result  says  that  we  can 
find  optimal  solutions  for  problems  (P)  and  (f>)  in  O(n’)  arithmetic  operations  once  the 
duality  gap  at  a  point  wk  generated  by  algorithm  3.1  becomes  sufficiently  small. 

Proposition  4.6:  Let  w  =  (a,  y,  c)  be  a  point  m  the  ser  W  satisfying 


x < 


( m  +  n)~ 


(4.25) 


Then  we  can  find  a  point  w*  —  (a*,  y* ,  :*)  in  no  more  than  O (n3)  arithmetic  operations, 
such  that  a*  and  vc*  =  (.v*,  y  *,  r*)  solve  problems  (P )  and  ( D )  respectively. 

This  result  can  be  proved  by  slightly  modifying  the  arguments  of  lemma  2  of  [I], 
Using  this  result,  we  obtain 

Corollary  4.7:  If  the  initial  penalty  parameter  p°  satisfies  log  pn  =  O(L)  then  algorithm 
3  1  solves  problem  (P)  in  at  most  O(^nL)  iterations 

Proof:  Using  the  previous  proposition,  we  tan  set  r  =  2~2L/(m  +  n):  as  the  tolerance  for 
the  duality  gap  in  algorithm  3.1.  From  proposition  4.5,  we  immediately  conclude  the  vali¬ 
dity  of  this  corollary.  G 

In  section  5,  we  will  sec  that  the  initial  penalty  parameter  can  be  chosen  to  satisfy 
log  =  O(L).  One  possible  choice  for  the  approximation  s  =  (a,  :  )  on  step  2  of  the 
algorithm  3.1  is  to  use  exact  data,  that  is,  to  set  s.  on  the  k'h  iteration,  equal  to  j*.  With 


-  * » *  ■  ^  ‘j*  *  t’- .  ■ .  »■-  .*  .**.'•.*  ■« 

jaA jlA  jl/l  i  aAiJ.  A  v  -  -  V—  -  -  V—  VL  » 
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this  choice  o(  s  ,  we  have  the  following  resuli: 

Corollary  4.8:  Algorithm  3  I  solves  problem  (Pi  in  no  more  than  O(n^^L)  iterations. 

Proof:  At  every  iteration,  the  computational  elfort  is  majorized  by  the  calculation  of  the 
inverse  of  the  matrix  [A(Zl  +  A'  1  Q  f  'A 1 A  T\  which  requires  0(ny)  arithmetic  operations. 

By  corollary  4  7,  algorithm  3.1  terminates  in  at  most  O(VnZ-)  iterations.  These  two  obser¬ 
vations  immediately  concludes  the  proof  of  the  corollary.  □ 

In  the  next  subsection,  we  present  an  alternative  choice  for  the  approximation  that 
makes  possible  to  reduce  the  complexity  of  algorithm  3.1  to  0(n*L)  arithmetic  operations. 

4.1.  A  Good  Choice  for  .v  and  : 

In  this  subsection,  we  show  that  the  complexity  of  algorithm  3.1  can  be  reduced  to 
O(mL)  arithmetic  operations.  The  arguments  in  this  subsection  arc  the  same  as  for  the 
linear  programming  case  presented  in  section  5  of  j  12],  We  should  point  out  that  this  idea 
for  reduction  of  the  complexity  was  first  presented  in  Karmarkar  [7]  and  subsequently  in 
Gonzaga  [6]  and  Vaidya  [14],  The  reduction  basically  consists  of  using  a  direction  that 
approximates  the  "exact"  direction  calculated  from  using  "exact"  data,  that  is,  the  current 
iterate.  In  our  case,  an  approximate  direction  is  implicit  in  the  choice  of  the  approximation 
s  .  In  this  section,  we  show  that  by  choosing  the  approximation  s  conveniently,  a  reduction 
in  the  average  work  per  iteration  is  obtained.  The  choice  of  the  approximation  s  is  made 
by  an  updating  scheme  as  follows  (In  the  procedure  below,  k  stands  for  the  iteration  count): 

Updating  scheme  4.9: 

For  k  :=  0,  set  x  :=  x“  and  z  :=  z" 

For  k  >  0  do 

For  i  =  l . n  do 


//  one  of  the  following  holds 


I  a*  -  x,  I 

(a) - -  >  r 

I  x.  I 


„*  _  r 


(b)  — - - I  >  y 


then  set  x,  :  =  a*  and  :=  r* 


Some  observaiions  arc  in  order  at  this  point.  In  order  to  calculate  the  directions  Ax 


and  Ay  determined  by  system  (3.1).  we  need  to  calculate  the  inverse  of  the  matrix 


A(Z  +  XQ)~'XA 7  =  A(A"'Z  +  Q)‘  A 


(4.26) 


where  s  s  (*,  7)  represents  the  approximation  for  the  current  iteration.  Let  s*  and  Bk 


denote  the  approximation  s  and  the  matrix  given  by  (4.26)  respectively  at  the  k'h  iteration 


of  the  algorithm  3.1.  Also  let  Dk  denote  the  matrix  (Z*)  ]Xt.  We  show  next  that  if  the 


matrix  Dk  differs  from  the  matrix  Dk_{  by  exactly  /  diagonal  elements  then  the  computation 


of  Bk  1  can  be  carried  out  in  0(n~l )  arithmetic  operations  by  means  of  /  rank-one  updates. 


Let  £  =  Q  +  Dk_\  and  F  =  Dk  -  Dk_ ,.  Then  we  obtain 


Bt =  AE~'At 


(4.27) 


Bk  =  A  (E  +  F  )~'at 


(4.28) 


Obviously,  £  is  a  diagonal  matrix.  Denote  the  ilh  diagonal  element  of  the  matrix  F  by  /,. 


By  assumption,  exactly  /  diagonal  elements  /,  are  non-zero.  For  simplicity  of  notation,  we 


assume  that  that  these  elements  are  the  first  /  diagonal  entries  of  the  matrix  F.  Then  the 


matrix  F  can  be  written  as 


F  =  t'f.u'iu-)7 


A 


'  ‘  MX: 


where  u  denotes  the  n-vector  where  all  components  are  zero  except  the  /'*  component 
which  equals  one.  Let  £  be  defined  as 

E0=E 


Ej  =  +  fju’iu')1  .  j  =  1 .  I 


(4.29) 


Note  that  £,  =  £  +  F.  We  observe  that  the  matrices  E0 ,  £, . £,  are  positive  definite, 

and  hence  invertible  matrices.  Applying  the  well  known  Shermann-Morcison  formula  of 
Linear  Algebra  to  the  matrix  £  as  given  in  expression  (4.29),  we  obtain  for  ;  =  1 . / 


r-  “I  _  f  -  I 
L1  ~  LJ~  I 


\  ~  e;\u’(uWe;\ 

1  +  fJu’)TEj\u’ 


Recursively,  we  can  obtain  £,  as  follows. 


£/■'  =  -  Xs,v(vf 


(4.30) 


where  the  scalars  g,  and  the  /(-vectors  v*  ,  /  =  1 . /  are  generated  by  the  following  itera¬ 

tive  procedure. 


Procedure  4.11:  Given  EJ 1  then. 


For  i  =  1  . 1  do 


1  +  f Ju1  )T  Ej\u’ 


v'  =  e;\u> 


EJ '  =  EJ_\  -  gjv'(v’)7 

Since  £,  =  £  +  £  and  using  expressions  (4.27),  (4.28)  and  (4.30),  we  obtain 


=  Bk_x  -  j^gj(A  v’HAy')t 
)  =  i 


(4.31) 


We  can  also  use  the  same  process  described  above  to  find  the  inverse  of  the  matrix  Bk 


using  expression  (4.31)  and  the  matrix  already  calculated  in  the  prev  as  iteration  of 

the  algorithm.  We  note  that  the  procedure  above  involves  0(n2l)  arithmetic  operations. 

Next  we  provide  an  upper  bound  on  the  number  of  diagonal  element  changes  that 
occurs  on  the  matrix  Z~'X  during  K  steps  of  algorithm  3.1.  Note  that  the  i'h  diagonal  ele¬ 
ment  of  the  matrix  (Z  “'x  )  changes  only  when  inequality  (a)  or  (b)  of  the  updating  scheme 
4.9  is  satisfied. 

The  following  result  can  be  proved  by  using  the  arguments  in  section  5  of  [6]. 

Proposition  4.12:  Let  ( \* )f=0  be  a  sequence  of  n-vcctors  with  all  components  positive  and 
satisfying 


||vt  +  l  -  <  p  ,  *  =  0,  1,  . ..  ,K  -  I 

where  p  is  a  positive  constant  less  than  one.  Define  the  sequence  (v*)f=o  recursively  as 
follows.  Set  v°  :=  and  for  k  >  1  and  i  =  1 . n  let 


if 


I  v*  -  v*-'  I 


—  >  r 


v\ I 


v*-‘ 


otherwise 


where  y  is  a  positive  constant  less  than  one.  Let  V ,  be  the  set  of  indices  k  defined  as 


VK  = 


lv/  - 


lv/-'l 


>  y  .  1  <  k  <  K 


and  let  I  V,*  I  denote  its  cardinality,  that  is,  the  number  of  times  the  i'h  component  of  the 
sequence  (v*)*_0  changes.  Then  the  following  inequality  holds 


(1  -  p)  log(l  -  y) 


1  =  1 


As  a  consequence  of  this  resuli,  we  have  the  following  corollary. 

■  K  t  1C 

Corollary  4.13:  Let  ( x  )*_<,  and  (:  )L=0  be  the  sequences  generated  by  algorithm  3.1  and 
consider  the  approximation  s  =  (!,?)  defined  as  in  the  updating  scheme  4.9  Consider  the 
follow  ing  two  sets: 

(  \ 


SK  = 


I  X*  -  X,  I 
I  X.  I 


>  y 


V 

Then  the  follow  ing  inequalities  hold. 


\  <  k  <  K 

, 

' 

1  <  *  <  A' 


X  l-S^I  <  4.5s~nK 
<  =  i 

£  I  7"*l  <A.5<nK 
i  =  i 


Proof:  This  result  follows  immediately  by  using  relations  (4.21),  (4.22)  and  proposition 
4.12.  C 

Thus,  the  total  number  of  rank-one  updates  that  occurs  during  K  steps  of  algorithm  3.1 
is  on  the  order  of  O('fnK).  As  a  consequence  of  this  result,  we  have  : 

Corollary  4.14:  Algorithm  3.1  coupled  with  the  updating  scheme  4.9  solves  problem  (P)  in 
no  more  than  0(nyL)  arithmetic  operations 
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Proof:  From  corollary  4.7,  we  know  that  algorithm  3.1  finds  an  optimal  solution  of  problem 
(P)  in  0{\nL)  iterations  Corollary  4.13  implies  that  the  total  number  of  rank -one  updates 
is  then  of  the  order  of  0[nL).  Since  each  rank -one  update  involves  0(n2)  arithmetic 
operations,  the  total  number  of  arithmetic  operations  is  then  of  the  order  of  0(n2L).  This 
completes  the  proof  of  the  corollary.  C 

5.  Initialization  of  the  Algorithm 

In  this  section,  we  show  how  to  initialize  algorithm  3.1,  in  order  to  solve  any  convex 
Quadratic  Programming  problem.  The  approach  is  to  use  a  transformed  problem  equivalent 
to  original  one  that  satisfies  the  initial  condition  (3.3).  Therefore,  by  solving  the 
transformed  problem,  we  are  able  to  obtain  a  solution  for  the  original  problem.  We  also 
show  that  the  "size'  of  the  transformed  problem  is  of  the  same  order  of  the  ’’size"  of  the 
original  problem,  where  by  "size",  we  mean  the  constant  L  defined  in  relation  (4.23).  This 
fact  implies  that,  with  respect  to  the  data  of  the  original  problem,  the  complexity  obtained 
in  section  4  is  still  valid. 

Consider  the  convex  quadratic  programming  problem 

—  —7—  ]  ^  7  — 

(P  )  min  c  x  +  .»  Q  x 

s.t  A  x  =  b 

T  >  0 

I  where  A  is  an  mxn  matrix  which  has  full  row  rank,  Q  is  an  fixn  positive  semi-definite 

matrix  and  b  ,  c  are  vectors  of  length  m  and  «  respectively.  We  assume  that  the  entries  of 
the  vectors  b  and  c  and  the  entries  of  the  matrices  A  and  Q  are  integral.  Let 

I 

k 

* 

I 


L  =  log  ;  large:.;  absolute  \aluc  of  the  determinant 
o)  any  square  suhmutt ix  of  M 


+  log:  (  max  I  o  I  )  +  log:  (  max  I  h§  I  )  +  log^  (  m  +  ft  ) 
j  ‘ 


where  M  is  ihe  (m  +  h)  x  (rii-rh  j  matrix  delined  as 


Q  at 


The  Karush-Kuhn-Tucker  necessary  and  sullicicnt  condition  for  x  e  R"  to  be  a  solution  of 
(P  )  is  that  there  exists  a  vector  >  e  R'”  such  that  x  and  y  satisfy 


:=c+Qx~A  \>  0 


A  x  =  h 


x  >  0 


r  7I  =  0 


which  can  be  rewritten  as 


~Q  A  7 


0  J  {?) 


x  >  0  ,  r  £  0 


I7:  =  0 


(5. 2. a) 


(5.2.b) 


(5.2.c) 


A  well  known  result  from  Linear  Complementary  theory  is  that  if  the  system  (5.2)  has  some 
solution  then  it  has  a  solution  which  is  a  vertex  of  the  polyhedron  given  by  (5. 2. a)  and 
(5.2.b).  The  following  lemma  is  a  well  known  result  whose  proof  follows  from  an  immedi¬ 
ate  application  of  Cramer’s  rule. 

Lemma  5.1:  Let  w  =  (x,  y.T)  be  a  vertex  of  the  polyhedron  given  by  (5. 2. a)  and  (5.2. b). 
Then  the  coordinates  of  h  arc  rational  numbers  with  numerator  and  denominator  less  than 


-  24  - 


or  equal  to  1L 


We  observe  that  any  solution  (J,  y)  of  system  (5.1)  is  an  optimal  solution  of  the  dual 
of  problem  (P ),  which  is  the  problem  given  by 

( D )  max  -  *  vT Q  v  -f  bTx 
2 

—  Q  v  +  Ar y  <  c 

where  v  is  an  n  -vector  and  y  is  an  m-vector.  Let  n  =  n  +  2  and  m  =  m  +  1.  Let  A  =  2L 
and  A  >0  be  a  large  constant  which  will  be  specified  more  precisely  later.  Consider  the 
transformed  problem  as  follows. 

(P)  min  c7  x  +  ^  xT  Q  x  +  K  xn 

s.t.  A  x  +  (b  -  A Ae)x„  =  b 
eTx  +  +  A.v„  =  nX 

x  2:  0  ,  T„_,  >0  ,  x„  >  0 

where  x  =  (.v( . x„_2)T  is  an  (n-2)-vector  and  and  x„  are  scalars.  The  dual  prob¬ 

lem  corresponding  to  problem  (  P  )  is  given  by 


(D)  max  -  \tQ  v  +  bTy  +  (nX)ym 


-  Q  v  +  A7y  +  e  ym  <  V 


ym  zo 


(b  -  A Aefy  +  ym<  K 


where  y  =  (y , . >„.|)T  is  an  (m-l)-vector,  ym  is  a  scalar  and  v  is  an  (n-2)-vector. 

These  problems  can  be  recast  in  the  notation  of  problems  (P)  and  (£>)  of  section  2  as  fol¬ 
lows.  Let  x  -  (xT ,  x„)T e  R\  y  =  (y r,  ym)T e  Rm  and  v  =  (vT ,  v„_,,  v„)Te  R" 


Define  b  e  Rm.  c  e  R 


and  A  e  Rmx"  as  follows. 


h } 

f~] 

r 

c 

A  0 

b  -XAe 

rJ 

C  = 

0 

A  = 

er  1 

A 

Let  Q  e  R"*"  denote  the  block  diagonal  matrix  as  follows. 

Q  =  diag(Q ,  0,  0)  (5.4) 

With  these  notations,  we  can  then  rewrite  problems  ( P )  and  ( D )  as  in  section  2.  We  refer 
to  these  two  formats  interchangeably. 

In  the  following,  we  adopt  the  convention  to  denote  the  optimal  value  of  a  quadratic 
programming  problem  ( P )  as  val(P)  and  the  value  of  the  objective  function  of  (P)  at  a 
feasible  point  x  as  valp(x).  We  now  present  the  relationship  between  the  optimal  solutions 
of  the  transformed  problem  (P),  and  its  dual  (/)),  with  the  optimal  solutions  of  the  original 
problem  (P),  and  its  dual  (D),  respectively.  Before  stating  the  relation,  we  make  the  fol¬ 
lowing  observation.  If  problem  (P )  has  an  optimal  solution  then,  by  the  observations 
preceding  lemma  5.1,  a  pair  (  7#,  y  „ )  must  exist  which  solves  the  system  (5.2)  and  which  is 
a  vertex  of  the  polyhedron  defined  by  relations  (5. 2. a)  and  (5.2.b).  The  pair  (1^,  y„)  is 
considered  in  the  following  result. 

Lemma  5.2:  Assume  that  problem  (P)  has  an  optimal  solution  and  let  (T^,  y  )  be  as 
above.  Assume  that  the  cost  coefficient  K  in  the  cost  vector  of  problem  (P)  satisfies 

K>(b-kAe)Ty „  (5.5) 

Then  the  following  statements  hold 

(1)  The  common  optimal  value  of  problems  (P)  and  (£> )  is  equal  to  the  common  optimal 
value  of  problems  (P  )  and  (D). 

(2)  If  x  -  (jr, . x„)T  and  (v.  y)  =  «v, .  v„)r,  ( v, . ym)T)  ore  optimal  solutions  of 

problems  (P)  and  ( D )  respectively.  then  x„  =  0  and  ym  =  0.  Moreover. 


& 


ll» 


x  -  (a,,  ...  ,  x„_2)  and  (v,  y)  -  ((v,,  ...  ,  v„_2)r,  (y,,  ...  ,  ym_j)r)  are  optimal  solutions 
of  (P  )  and  ( D ). 

Proof:  We  first  prove  (1).  Consider  the  vectors  a„  e  R"  and  y*  e  R"  defined  as  follows. 


**  =  (* *  •  nk-eTx\ ,  0) 


=  (y*.o> 

Expression  (5.5)  and  lemma  5.1  imply  that  a^  and  (xt,  y\)  are  feasible  solutions  of  prob¬ 
lems  ( P )  and  ( D )  respectively.  Thus,  we  have 


val(P )  =  valp(x  *)  =  valp(x^)  >  val(P) 


val(D)  =  valD(xt.  yj  =  valD(x y*)  <  val(D) 


Since  val(D)  =  val(P)  and  val(D)  =  val(P),  relations  (5.6)  and  (5.7)  then  immediately 
imply  (1).  Moreover,  x *  and  (a#,  y^)  are  optimal  solutions  for  (P)  and  (D)  respectively. 

We  now  prove  (2).  Since  x*  and  (v,  y)  form  a  pair  of  primal  and  dual  optimal  solu¬ 
tions  for  problems  (P)  and  ( D )  respectively,  they  must  satisfy  the  complementary  slackness 
condition  (c.f.  proposition  2.2  of  section  2).  In  particular,  we  have 

(->«)  (nk-eTxJ  =  0  (5.8) 

But  lemma  5.1  implies  that  eTx *  <  n k  =  (n-2)k  <  nk.  Therefore,  (5.8)  implies  that 
ym  =  0.  But  this  implies  that  (v,  y)  is  feasible  to  (D)  and  that  va/D(v,  y)  =  val^(v,  y). 
Statement  (1)  above  then  implies  that  (v,  y)  is  optimal  for  ( D ). 

Arguing  with  the  pair  a  and  (a+,  y*)  in  a  similar  way,  we  conclude  that  x„  =  0  and 
that  .7  is  optimal  for  (P ).  This  proves  (2).  □ 


Lemma  5.3:  Assume  that  problem  ( P )  has  an  optimal  solution  and  let  (x  y*)  be  as  in 
lemma  5.2.  Then  K  =  2^  satisfies  relation  (5.5). 

Proof:  This  lemma  follows  straightforwardly  from  the  definition  of  L  and  from  lemma  5.1. 


In  view  of  lemma  5.3.  from  now  on,  we  lei  K  =  2?1 .  Consider  the  constant  L  defined 
as  in  relation  (4.23)  and  (4.24).  The  following  observations  are  easily  shown. 

(1)  From  the  definition  of  Q,  A,  b  and  c  given  by  expressions  (5.3)  and  (5.4),  it  immedi¬ 
ately  follows  that 


L  <  L 


(2)  The  largest  absolute  value  of  the  determinant  of  any  square  submatrix  of  A  is  at  most 
(in  +  h)  23*\ 

(3)  max  1  b,  I  <  nl^  and  max  1  c  I  £  22^ . 

i  =  l . m  )=\ . n 

(4)  Statements  (2)  and  (3)  implies  that  L  <  9L . 

We  now  verify  that  problem  (P)  satisfies  assumption  2.1  of  section  2.  Assumption  (c) 
is  obviously  satisfied  since  A  was  assumed  to  have  full  row  rank.  We  verily  assumptions 
(a)  and  (b)  jointly  by  exhibiting  a  point  w °  =  (  v",  y°,  z°)  which  is  in  the  set  W  defined  in 

section  2  and  satisfying  the  criterion  of  closeness  (3.3).  Let  x°  s  (A . A,  l)re  R".  Let 

Qr  j  =  1 . n  denote  the  j'h  row  of  the  matrix  Q  and  let  y°  =  (0,  ...  ,  0,  -p°lk)T  e  Rm 

where  p°  satisfies 


.  £(Ac,  +  A Q/x0)2  +  (c„  +  Q„tx")2 
j  =  i 


Let  r°e  R"  denote  the  slack  vector  c  +  Qx"  -  A  Ty°  for  the  dual  (D)  corresponding  to  the 
pair  v  =  x°  and  y  -  y°.  Since  Ary°  =  p°lk  (1 . l,A)r  e  R\  it  is  easy  to  verify  that 


5>>V  ■  »0) =  X(Af; +  xQ,Tx"r  +  +  e/*0)2 

;=1  y=l 

and  hence  that  the  criterion  of  closeness  (3.3)  is  satisfied  due  to  expression  (5.9).  Since 
logA  =  L  <  L,  it  is  straightforward  to  verify  that  the  logarithm  of  right  hand  side  of  (5.9) 
is  on  the  order  of  O(L).  Therefore  p°  can  be  chosen  to  satisfy  log  p°  =  O(L).  By  the 
convergence  results  of  section  4  and  5,  it  follows  that  problems  ( P )  and  (D)  can  be  solved 
in  at  most  O(n^L)  arithmetic  operations.  Finally,  the  main  result  of  this  section  is  as  fol- 


Proposition  5.4:  Problem  (P )  can  be  solved  in  at  most  0(nL)  arithmetic  operations. 

Proof:  Applying  algorithm  3.1  to  problem  (P),  we  obtain  vectors  x  =  (xj,  ...  ,  xn)  and 

>  =  (>! . ym)  such  that  x  and  (x,  y)  are  optimal  solutions  for  problems  (P)  and  ( D ) 

respectively.  Consider  the  following  two  possible  cases. 

(i)  If  x„  =  0  and  ym  =  0  then  problem  (P),  and  consequently  (£>),  has  an  optimal  solution. 

Indeed,  if  we  let  x  =  (x, . x„_2)T  and  y  =  (y, then  x  and  ( x ,  y)  are  feasi¬ 

ble  solutions  for  ( P )  and  (D)  respectively.  Statement  (2)  of  lemma  5.2  then  implies  that  x 
and  (x,  y )  are  optimal  solutions  for  ( P  )  and  ( D )  respectively. 

(ii)  If  either  xn  *  0  or  ym  *  0  then  lemma  5.2  implies  that  (P)  is  either  unbounded  or 
infeasible.  In  this  case  we  solve  the  LP  problem  obtained  by  replacing  the  objective  func¬ 
tion  of  problem  (P)  by  the  linear  function  Kx„.  If  the  resulting  optimal  solution  of  this 
problem  satisfies  =  0  then  ( P )  is  unbounded.  Otherwise,  (P )  is  infeasible. 

By  corollary  4.14,  the  computation  above  can  be  carried  out  in  at  most  CHn^L)  arithmetic 
operations.  Since  n  =  n  -  2  and  L  <  9L ,  the  total  number  of  arithmetic  operations  is  on 
the  order  of  0(H3L ).  □ 


6.  Remarks 


The  following  observations  are  in  order: 

(1)  The  purpose  of  this  paper  is  to  present  a  theoretical  result.  Thus  in  order  to  simplify  the 

presentation,  we  constructed  (i  =  p(l  —  5  /Vn).  Obviously,  one  can  use  /}  which  is  less 
than  or  equal  than  the  above  one,  but  still  satisfying  (b)  of  theorem  4.3  and  relations  (4.21) 
and  (4.22).  In  this  way,  one  can  accelerate  the  convergence  of  the  algorithm. 

(2)  Additional  improvements  in  actual  implementation,  which  are  possible,  such  as  more 
judicious  selection  of  6,  S  and  y,  together  with  actual  test  results,  are  the  subject  of  a  forth¬ 
coming  paper. 

(3)  With  the  necessary  modifications,  the  results  of  this  paper  are  also  valid  in  the  case  that 
the  matrix  Q  is  only  assumed  to  be  positive  semi-definite  on  the  affine  space 
{  x  I  Ax  =  b  ).  All  the  duality  results  remain  true  if  we  add  the  constraint  Ax  =  b  to  the 
dual  problem  (D).  The  formulas  for  the  directions  Ax  and  Ay  given  in  section  3  do  not 
necessarily  hold  in  this  case.  We  leave  it  to  the  reader  to  carry  out  the  necessary 
modifications.  Finally,  we  note  that  the  convergence  results  of  section  4  follow  without  any 
modifications. 

(4)  It  is  well  known  that  a  Linear  Complementarity  problem  with  positive  semi-definite 
matrix  can  be  reduced  to  an  equivalent  convex  Quadratic  Programming  problem  and  vice- 
versa  (  c.f.  (11]  ).  Thus,  the  algorithm  presented  in  this  paper  can  be  used  to  solve  Linear 
Complementarity  problems  with  positive  semi-definite  matrices.  At  the  lime  of  writing  this 
paper,  we  were  informed  of  a  recent  paper  by  Kojima  ct.  al.  [10]  which  present  an  algo¬ 
rithm  for  solving  Linear  Complementarity  problems  with  positive  semi-definite  matrices. 
They  obtained  the  same  complexity  as  the  one  achieved  in  this  paper. 
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